Home Range Analysis
Advanced Data Analytics

library("adehabitatHR")
library("data.table")
library("ggfortify")
library("grid")
library("move")
library("moveVis")
library("pbapply")
library("plotly")
library("rgdal")
library("sp")
library("tidyverse")
library("viridis")
library('ggmap')
library("OpenStreetMap")
library("gpclib")
library("maptools")

1 Dataset

This data comes from my thesis research from this past spring/summer. This dataset contains only my study site on the street behind Dr. Schillers’ property, which contains GPS points for 3 individual male Louisiana Waterthrushes.

Louisiana Waterthrush

Importing dataset with read.csv:

data <- read.csv("./Data/lowaSCH.csv")

qaqc_plot <- ggplot() + geom_point(data=data, 
                                   aes(location.long,location.lat,
                                       color=individual.local.identifier)) +
                        labs(x="Longitude", y="Latitude") +
                        guides(color=guide_legend("Individual"))

ggplotly(qaqc_plot)

With plotly, similar to leaflet, we have the ability to examine the spread of the data points and additional information from various columns. From this plot we can see that there are three individuals; m1, m2, and m3. Already you can see how their territories follow pretty tightly along the stream and tributary at this site.

While we could continue with the current dataset, any analysis would calculate home range for the entire population rather than the individual.

Creating a function and using the lapply command to apply a function over a list or vector dataset. Specifically, this function will take the original dataset, split it into separate files based on the individual identifier, and create new *.csv files using the identifier as the filename.

lapply(split(data, data$individual.local.identifier), 
       function(x)write.csv(x, file = paste(x$individual.local.identifier[1],".csv", sep = ""), row.names = FALSE))

files <- list.files(path = ".", pattern = "[m]+[0-9]+", full.names = TRUE)

In the list.files command above, path = "." informs the locations, in this case the root directory, pattern = describes the way the files are named, in this case “m” (for male) followed by a number between 0-9, and full.names describes how the files will be listed.

2 Analysis

2.1 Imagery

Although my data contained only longitude/latitude values, I converted them to UTM and back again to longitude/latitude coordinates to practice using the code if I ever receive data that only has UTM values.

utm_points <- cbind(data$utm.easting, data$utm.northing)
utm_locations <- SpatialPoints(utm_points, 
                 proj4string=CRS("+proj=utm +zone=16 +datum=WGS84"))
proj_lat.lon <- as.data.frame(spTransform(
                utm_locations, CRS("+proj=longlat +datum=WGS84")))
colnames(proj_lat.lon) <- c("x","y")
raster <- openmap(c(max(proj_lat.lon$y)+0.001, min(proj_lat.lon$x)-0.001), 
                  c(min(proj_lat.lon$y)-0.001, max(proj_lat.lon$x)+0.001), 
                  type = "bing")
raster_utm <- openproj(raster, 
              projection = "+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs")

In the script above, utm_point is an x,y derived from the primary dataset, utm_locations set the projection to UTM Zone 16, proj_lat.lon converted the UTM points to longitude/latitude, raster uses the min/max x,y data to create a bounding box to retrieve the aerial imagery, and raster_utm reprojected the imagery back to UTM Zone 16 for Tennessee. Next, autoplot.OpenStreetMap will display the raster image file with the UTM locations as an overlay.

autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() +
  theme(legend.position="bottom") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_point(data=data, aes(utm.easting,utm.northing,
             color=individual.local.identifier), size = 3, alpha = 0.8) +
  theme(axis.title = element_text(face="bold")) + labs(x="Easting",
        y="Northing") + guides(color=guide_legend("Identifier"))

3 Home Range Analysis

In the section above we use the lapply command to loop a function used to separate the original dataset into individual files. This is a useful tool, however, when the function loops through dozens or even hundreds of files, the process can take a long period of time to complete. Using the pblapply command adds a progress bar (i.e. pb) to the process which provides an estimated time for completion of the function. We will use a similar process to the one above, using the pblapply command to run MCP analysis on the individual *.csv files.

3.1 Minimum Convex Polygon

A description of the steps within the following code will be discussed following the output. The process works by establishing the function and then running pblapply referring back to the files we created in the dataset portion of this exercise.

I had some issues getting “pblapply” to run, so I had to install the package “gpclib” and specify type=“source” in order to bypass whatever error I was getting.

I also changed the point and polygon colors in order to see them better on the dark background (especially since I couldn’t get it to zoom in any closer to the points).

if (!require(gpclib)) install.packages("gpclib", type="source")
gpclibPermit()

[1] TRUE

mcp_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$utm.easting)
  y <- as.data.frame(data$utm.northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  mcp.out <- mcp(xy, percent=100, unout="ha")
  mcp.points <- cbind((data.frame(xy)),data$individual.local.identifier)
  colnames(mcp.points) <- c("x","y", "identifier")
  mcp.poly <- fortify(mcp.out, region = "id")
  units <- grid.text(paste(round(mcp.out@data$area,2),"ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  mcp.plot <- autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=mcp.poly, aes(x=mcp.poly$long, y=mcp.poly$lat), color="#ff5e00",fill="#c34800", alpha=0.7) +
    geom_point(data=mcp.points, aes(x=x, y=y), color = "#ff9b00", size=3) + 
    labs(x="Easting (m)", y="Northing (m)", title=mcp.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  mcp.plot
}

pblapply(files, mcp_raster)

  |                                                         | 0 % ~calculating  
  |===================                                      | 33% ~01s          
  |======================================                   | 67% ~00s          
  |=========================================================| 100% elapsed=01s  
[[1]]

[[2]]

[[3]]

3.2 Kernel-Density Estimation

I kept getting errors with this map, saying that “The grid is too small to allow the estimation of home-range. You should rerun kernelUD with a larger extent parameter”.

After trying all kinds of code fixes, what finally got this to run was changing the “getverticeshr(kde, 95)” from 95 to 30. I’m not sure why this is the magic number that worked, but it did give me some home range estimates as you can see in the maps above.

kde_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$utm.easting)
  y <- as.data.frame(data$utm.northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  kde<-kernelUD(xy, h="href", kern="bivnorm", grid=100)
  ver <- getverticeshr(kde, 30)
  kde.points <- cbind((data.frame(data.proj@coords)),data$individual.local.identifier)
  colnames(kde.points) <- c("x","y","identifier")
  kde.poly <- fortify(ver, region = "id")
  units <- grid.text(paste(round(ver$area,2)," ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  kde.plot <- autoplot.OpenStreetMap(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=kde.poly, aes(x=kde.poly$long, y=kde.poly$lat), color="#0d9be8", fill="#0d9be8", alpha = 0.8) +
    geom_point(data=kde.points, aes(x=x, y=y), color="#6ee1f0", size=4) +
    labs(x="Easting (m)", y="Northing (m)", title=kde.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  kde.plot
}

pblapply(files, kde_raster)

  |                                                         | 0 % ~calculating  
  |===================                                      | 33% ~01s          
  |======================================                   | 67% ~00s          
  |=========================================================| 100% elapsed=01s  
[[1]]

[[2]]

[[3]]

3.3 Brownian Bridge Movement

Although this appears to be working backwards, in the previous two examples you have seen how to create a function that can be looped to run the analysis for a list fo files. For this analysis, we will create a conventional script for a single individual. Portions of the script that were not previously described will be discussed following the analysis.

LOWAsch <- read.csv("./Data/lowaSCH.csv")

date <- as.POSIXct(strptime(as.character(LOWAsch$timestamp),"%Y-%m-%d %H:%M:%S", tz="US"))
LOWAsch$date <- date
LOWAsch.reloc <- cbind.data.frame(LOWAsch$utm.easting, LOWAsch$utm.northing,
                                as.vector(LOWAsch$individual.local.identifier),
                                as.POSIXct(date))
colnames(LOWAsch.reloc) <- c("x","y","id","date")
trajectory <- as.ltraj(LOWAsch.reloc, date=date, id="LOWAsch")

Error in as.ltraj(LOWAsch.reloc, date = date, id = "LOWAsch") : 
  non unique dates for a given burst

Issue w/ dates

3.4 Animate Trajectory Data

The final analysis we will perform in this exercise is a visual animation of the relocation data. Using the data from above you can plot(trajectory) and see a plot of the relocation information for an individual. However, we can use the move and moveVis packages to create an animation of the relocations.

We need to begin by creating a move object containing relocations (x,y), time and date information (time), a projection string, individual identifier (animal), and sensor type (sensor). See the description for move() in the help menu for optional information.

lowa_move <- move(x=LOWAsch$location.long, 
#lowa2_move<-move(df$location.long, df$location.lat,
             y=LOWAsch$location.lat, 
             time=as.POSIXct(LOWAsch$timestamp, 
                             format="%m/%d/%Y %H:%M",tz="UTC"), 
             #time=as.POSIXct(df$timestamp,format="%m/%d/%Y %H:%M",tz="UTC"),
               proj=CRS("+init=epsg:32615"))

movement <- align_move(lowa_move, res = "max", digit = 0, unit = "secs")

With the data now on a uniform time scale we can create the frames and animation for the relocations. For this step I will use a basemap from MapBox which require token access through the use of their API. To do this you need to register with MapBox and create an access token. Then create a .Renviron file in your project folder. Copy the token information from MapBox and create an object in the .Renviron file such as map_token = 'paste token here' and add map_token = Sys.getenv('map_token') to the script below. However using the get_maptypes() script you can see there are various map services and map types that can be used. A simple output would be to use map_service = 'osm' (OpenStreetMaps) and map_type = 'topographic' or other map types available by viewing get_maptypes('owm'). when using a basemap without token access the map_token option can be removed from the script below.

frames <- frames_spatial(movement, path_colors = "red",
                         map_service = "osm", 
                         map_type = "topographic",
                         alpha = 0.5) %>% 
  add_labels(x = "Longitude", y = "Latitude") %>%
  add_northarrow() %>% 
  add_timestamps(movement, type = "label") %>%
  add_scalebar(distance = 2) %>% 
  add_progress()

Checking temporal alignment...
Processing movement data...
Approximated animation duration: ≈ 1.12s at 25 fps for 28 frames
Retrieving and compositing basemap imagery...

  |                                                         | 0 % ~calculating  

Error in curl_download(url = url, destfile = file) : HTTP error 404.

animate_frames(frames, fps = 5, overwrite = TRUE,
               out_file = "./moveVis-2021lowa.gif")

Error in animate_frames(frames, fps = 5, overwrite = TRUE, out_file = "./moveVis-2021lowa.gif") : 
  object 'frames' not found